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Puiseux series are power series in which the exponents can be fractional and/or 
negative rational numbers. Several computer algebra systems have one or more 
built-in or loadable functions for computing truncated Puiseux series ~ perhaps 
generalized to allow coefficients containing functions of the series variable that are 
dominated by any power of that variable, such as logarithms and nested logarithms 
of the series variable. Some computer-algebra systems also offer functions that 
Q ■ can compute more-general truncated recursive hierarchical series. However, for all 

of these kinds of truncated series there are important implementation details that 
haven't been addressed before in the published literature and in current implemen- 
K^ ' tations. 

For implementers this article contains ideas for designing more convenient, cor- 
^f) [ rect, and efficient implementations or improving existing ones. For users, this article 

is a warning about some of these limitations. Many of the ideas in this article have 
Cn . been implemented in the computer-algebra within the TI-Nspire calculator, Win- 

Z^ I dows and Macintosh products. 

, , 1 Introduction 

X 

^ . 

C^ . Here is a conversation recently overheard at a car-rental desk: 

Customer: "I followed your directions of three right turns to get on the high- 
way, but that put me in a fenced corner from which I could only turn right, 
bringing me back to where I started!" 

Agent: "Make your first right turn after exiting the rental car lot." 

The original directions were correct, but incomplete. 

The same was true of published algorithms for truncated Puiseux series. After read- 
ing all that I could find about such algorithms, I implemented them for the computer 
algebra embedded in the TI Nspire™ handheld graphing calculator, which also runs on 
PC and Macintosh computers. Testing revealed some incorrect results due to ignorance 
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about some important issues. Results for other series implementations, reveal that their 
implementers have made similar oversights. 

It required a significant effort to determine how to overcome these difficulties. This 
article is intended as a warning for users of implementations that exhibit the flaws - 
and as suggestions to implementers for repairing those flaws or avoiding them in new 
implementations. Many of the ideas here are implemented in TI-Nspire. 

Additional issues for truncated and infinite series are described in a sequel to this 
article tentatively titled "Series crimes" |13) . 

For real-world problems, exact closed-form symbolic solutions are less frequently ob- 
tainable than are various symbolic series solutions. Therefore in practice, symbolic series 
are among the most important features of computer algebra systems. 

Almost all computer algebra systems have a function that produces at least trun- 
cated Taylor series. Iterated differentiation followed by substitution of the expansion 
point provides a very compact implementation. However, it can consume time and data 
memory that grows painfully with the requested order of the result. Knuth [5j presents 
algorithms that are significantly more efficient for addition, multiplication, raising to a 
numeric power, exponentials, logarithms, composition and reversion. He also suggests 
how to derive analogous algorithms for any function that satisfies a linear differential 
equation. Silver and Sullivan [lOj additionally give such algorithms for sinusoids and hy- 
perbolic functions. Brent and Kung [1] pioneered algorithms that are faster when many 
non-zero terms are needed. 

If for a truncated Puiseux series expanded about z = the degree of the lowest-degree 
non-zero term is a and g is the greatest common divisor of the increments between the 
exponents of non-zero terms, then the mapping z^ — )■ t(^~°)/5 can be used, with care, to 
adapt many of the Taylor series algorithms for Puiseux series. 

As described by Zipple [HI [15], many Puiseux series implementations have generalized 
them to allow coefficients that contain appropriate logarithms and nested logarithms that 
depend on the series variable. Geddes and Gonnet [3j, Gruntz j3], and Richardson et. 
al. [9\ give algorithms for more general truncated hierarchical series that also correctly 
prioritize essential singularities and perhaps-nested logarithms in coefficients. Koepf [7j 
implemented inhnite Puiseux series in which the result is expressed as a symbolic sum 
of terms: The general term in the summand typically depends on the summation index 
and powers of the series variable. Some implementations compute Dirichlet, Fourier, or 
Poisson series. Most of the issues described in this article are relevant to most kinds of 
truncated series, and some of the issues are also relevant to infinite series. 

To obtain a series for / (w) expanded about w = Wq with wq finite and non-zero, we 
can substitute w -^ z + wq into / (w) giving g {z) , then determine the series expansion 
of g (z) about z = 0, then back substitute z -^ w — Wq into that result. However, if 
the series is to be used only for real w < Wq then using instead w -^ Wq — z might give 
a result that more candidly avoids unnecessary appearances of i - particularly if / (w) 
contains logarithms or fractional powers. 

To obtain a series for / (w) expanded from the complex circle of radius oo, we can 
substitute ( — )■ 1/z into / (w) giving h (z), then determine the series expansion of h (z) 
about z = 0, then back substitute z — )■ 1/w into that result. 



A proper subset of the complex circle at infinity can be expressed by an appropriate 
constraint on the series variable, such as 

series (/ (w) , w = — oo, . . .) — )■ H { — 

\w 

where 

Hi.) = senes(/(i),. = 0,...)|.<0. 

Therefore without loss of generality, the expansion point is 2 = throughout the remain- 
der of this article with, z = x + iy = re*^ where r > 0, — tt < 6 < tt, and x,y ^W. 

Also, wherever braced case constructs occur, the tests are presumed to be done using 
short-circuit evaluation from top to bottom to avoid the clutter of making the tests 
mutually exclusive. 

2 The disorder of order 

"What we imagine is order is merely the prevaihng form of chaos." 

- Kerry Thornley 

Truncated series functions usually have a parameter by which the user requests a 
certain numeric "order" for the result. Existing implementations treat this request in 
different ways, some of which are significantly more useful than others. 

2.1 Render onto users what they request 

"Good order is the foundation of aU things." 

- Edmund Burke 

"Mathematics is the art of giving the same name to different things." 

- Jules Henri Poincare 

Unfortunately, the word "order" is used in too many ways in mathematics. Relevant 
definitions used in this article are: 

Definition. The exact error order of a truncated series result expanded about z = is 
T if the error is O {z'^) but the error isn't {z'^). The exact error order of an exact series 
result expanded about z = is +00. 

Remark. Knuth [6j introduced the convenient notation O {z'^) to denote exact-order r in 
z. In comparison to O {z'^), [z^) avoids discarding valuable information when we also 
know that a result isn't {z'^). 

Definition. The degree of a truncated Puiseux series with respect to z expanded about 
z = is the largest exponent of z that occurs outside of any argument of any O (. . .), 
o (. . .), or (. . .) term that is included in the result. 



It is unreasonable to request a degree because, for example, there is no way for 

series (cos 2;, z = 0, degree = 1) 

to return a result of degree 1. It is also unreasonable to request an exact error order 
because, for example, series [cosz, z = 0, O {z^)) can't return a series having error O (z^). 

Definition. If a series-function order-argument r denotes a request that the result be 
■ ■ ■ + {z'^), then the degree of an as-requested big-0 result should be the largest degree 
that satisfies 

degree < r < exact error order. (1) 

Remark. It seems likely that more often users prefer to specify the highest degree term 
they want to view rather than the lowest degree term they don't want to view. Thus 
most users would prefer that the series function parameter r denotes a request for a 
result that is ■ ■ ■ + o [z'^)- Therefore: 

Definition. If a series-function order-argument r denotes a request that the result be 
• ■ ■ + o {z'^), then the degree of an as-requested little-o result should be the largest degree 
that satisfies 

degree < r < exact error order. (2) 

The Maxima, Mathematica® and TI-Nspire™ truncated Puiseux series functions use 
this little-o interpretation of a numeric parameter r. For example, glossing over their 
input and output syntax differences, they all give 



sinz „ ^A 11 z^ z^ 

;z3 ' ' y z^ Q 120 5040 



series ( — ^, ^ = 0, 5 I -)■ —-^ + 77^-777777- (j) 



To this Maxima appends "+ ■ ■ ■" and Mathematica appends "+ 0[2;]^". 

It is easier to implement an interpretation in which the order parameter r denotes that 
the inner-most sub-expressions are computed to {z'^) and the final result is computed to 
whatever order that yields. For reasons described below, that can and often does lead to 
a result that is {z'^) with a n that is smaller or occasionally larger than r. For example, 
it would omit the last term of result ([3]). For such an implementation it is essential to 
display an error term because otherwise: 

1. If exact order < requested little-o order, then the result doesn't reveal that it is 
less accurate than requested, which can be disastrous. 

2. If degree > requested little-o order, then the user must notice and perhaps somehow 
truncate the excess terms to use the result in further calculations as intended. 

However, even if there is an error term indicating a result that doesn't have the requested 
accuracy, this design is inconvenient for users because: 

1. Users often don't notice the deficient or excessive accuracy. 

2. Users who notice excessive order, must perhaps somehow truncate the excess terms 
to use the result in further calculations as intended. 



3. Users who notice deficient order are forced to iteratively guess tlie order argument 
to use in series(. . .) to obtain sufficient accuracy - then perhaps somehow truncate 
a result that exceeds the desired order. 

4. If the user is another function, then that function should test the returned order 
and correct it if necessary by iterative adjustment and/or truncation. This is a 
requirement that might not occur to many authors of such functions - particularly 
those who aren't professional computer-algebra implementers. 

It is more considerate, reliable and efficient to build any necessary iterative adjustment 
and/or truncation into the series(. . .) function rather than to foist it on all function im- 
plementers and top-level users. It isn't prohibitively harder to implement an as-requested 
result. 

2.2 How to deliver as-requested order 

Definition. If an infinite Puiseux series is 0, then its dominant term is 0. Otherwise the 
dominant term is the lowest-degree non-zero term. 

Definition. If an infinite Puiseux series is 0, then its dominant exponent is oo. Otherwise 
the dominant exponent is the exponent of the dominant term. 

Remark. Some authors call the dominant exponent the valuation or valence, but other 
authors confusingly call it the order. 

A typical truncated-Puiseux-series implementation recursively computes series for the 
operands of each operator and the arguments of each function, combining those series 
according to various algorithms. Table [H lists the dominant exponent of a result and the 
operand orders that are necessary and sufficient to determine a result to o (2;*^). In that 
table a result dominant exponent of —00 signifies an essential singularity. 

As indicated there, cancellation of the dominant terms of series U and V can cause 
the dominant exponent of U ±V to exceed min {a, (5) when a = /3, such as for 

series (e^ — cos z, 2; = 0, o(z)) — )■ {{1 + z + o{z)) — {1 + o{z)) 

-)■ z + o{z). 

If the coefficient domain has zero-divisors, such as for modular arithmetic or fioating- 
point with underfiow, then the dominant exponent of UV can exceed a + (3, and the 
dominant exponent of U"' can exceed 7a. 

Unfortunately, most of the entries in column 3 require us to know the dominant 
exponents of the operands, perhaps also together with a dominant coefficient c and the 
exponent a of the next non-zero term, if any. Therefore we need this information before 
computing the operand series to the correct order, but we don't have this information 
until after we have computed the first term or two of the operand series. 

One way to overcome the difficulty is as follows: We can guess the dominant expo- 
nent of the operands by using a function written according to rewrite rules such as the 



following, which are heuristically motivated by the second column of Table [T] 

gViessDE{z,z) — ?■ 1, 

guessDE (m + V, z) — )■ min (guessDE (m, z) , guessDE (t>, z)) 

guessDE {uv, z) — )■ guessDE (m, z) + guessDE(t>, z), 

guessDE (m^, z) —> k guessDE {u, z) , 

guessDE (e", z) ^ 0, 

guessDE (sin u, z) — )■ max (0, guessDE (m, z)) , 

guessDE (m —1,2;), if m(0) = 1, 
0, otherwise, 



guessDE (In-u, z) 



{guessDE [u — i, z) , if m(0) = i, 

guessDE (m + i, z) , if m(0) = — i, 

max (0, guessDE (u, z)) , otherwise, 
guessDE (m, z) | u is independent of z — > 0. 

If using the guess results in computing more terms than necessary, then we should 
truncate the excess. If using the guess doesn't produce the required order but reveals the 
dominant term (and where needed the next non-zero term), then we know precisely the 
necessary and sufficient order to request for recomputing the series operands. 

If using the guess doesn't reveal this information, then when there is only one function 
argument we can iteratively increase the guess, starting with an initial increment 6 > 0. 
For each iteration we can double the increment added to the initial guess. This way, 
in a modest multiple of the time required for the last iteration, the process terminates 
successfully or by resource exhaustion. 

Resource exhaustion can be caused by an undetected essential singularity, insufficient 
simplification of an operand expression, or undetected constancy around the expansion 
point, such as for |x + l| + |a:; — 1| atx = 0, which can be more candidly expressed as 

-2x X < —1, 
1 -l<x<l, 

2x otherwise. 

To increase the likelihood of the first increment 5 being sufficient to expose the first 
non-zero term or two but not prohibitively more terms than needed, we can use a function 
that guesses the increment between the exponents of the first two non-zero terms. Let u 
and V be expressions with a = guessDE(u, z) and (3 = guessDE(v, z). Then the following 



ordered rewrite rules are examples for such a function: 



guessinc {z, z) — )■ 0, 
guessinc {v? , z) — )■ guessinc (m, -z) , 

guessinc (-u, z 
guessinc {uv, z) ^ I guessinc (f , z 



guessinc (In u) — )■ 

guessinc (e", z) — )■ 

guessinc (sin -u, z) — )> 

guessinc (cos u, z) — )■ 

guessinc (arctan-u, z) — )> 



guessinc {v,z) = 
guessinc (m, z) = 
min (guessinc (m, z) , guessinc (f , z)) , otherwise, 

guessinc {u — l,z) if -u (0) = 1, 
guessinc (m, z) otherwise, 

guessinc (m, z) if a = 0, 
\q.\ otherwise, 

2 I a I if guessinc (-u, z) = 

guessinc (-u, z) otherwise, 

guessinc (m, z) if a = 0, 
2 I a I otherwise, 

—a if (5 < 0, 

2(5; if q; > A guessinc (-u, 2;) = 

guessinc {u — m(0), z) if ■u(O) = i M ■u(O) = — i, 

^guessinc (-u, z) otherwise. 



// Same for sinh u 
II Same for cosh u 



guessinc (arctanh-u, z) — )> < 



guessinc (arccosh m, z) — )■ 



—a 
2a 

guessinc 
^ guessinc 

{2|a| 
guessinc 
guessinc 

{2|tt| 
guessinc 
guessinc 



a 

-2a 
guessinc 
^ guessinc 



if q; < 0, 

if q; > A guessinc (-u, 2;) = 0, 
m-m(0),2) if m(0) = 1 V m(0) = -1, 
■u, z) otherwise, 

if guessinc (-u, 2;) = 0, 
u,z)l2 ifM(O) = 1 V m(0) = -1, 
■u, z) otherwise, 

if guessinc (-u, 2;) = 0, 
u,z) 12 ifM(O) = i V m(0) = -i, 
■u, z) otherwise, 

if a > 0, 

if a < A guessinc (m, 2;) = 0, 
M,2)/2 ifM(O) = 1 V m(0) = -1, 
■u, 2;) otherwise. 



guessinc (m, 2;) | m is independent of 2; — )■ 0. 



The treatment of sums is more easily stated procedurally: Let s be a sum and a = 
guessDE(s, 2). Let S_ be the set of all terms in s that have a as the guess for their 
dominant exponent, and let S be the set of all the other terms. If S is empty, let 7 = 0. 
Otherwise let 7 be the minimum guessed dominant exponent of the terms in S. Let (5 = 
if is the guessed increment for all of the terms in S_. Otherwise let S be the minimum 
non-zero guessed increment in S_. Then the guessed increment for the sum is 



min f (5, 7 — cr 



Because guessDE(. . .) might return an incorrect guess for a dominant exponent, it 
is possible for this guessinc (. . .) to return when the increment is actually positive. 
Therefore we can instead guess an increment of 1 if a top-level invocation of guessinc (...) 
returns 0. 

The required argument order in Table [1] for the inverse trigonometric and inverse 
hyperbolic functions depends on whether the dominant exponent a of the argument 
u is positive, negative, or with the corresponding coefficient being a branch point. 
Evaluating u (0) can help us decide this: If u (0) is a branch point, then we can use the 
above iteration scheme on u — u (0) to obtain the required order o (2;'^"'"'^) . If u (0) = 0, 
then a > 0. If m (0) is otherwise finite and non-zero, then « = 0. Either way, the 
required order for -u is o (z'^) . Otherwise, either a negative dominant exponent or a 
coefficient having a logarithmic singularity caused u (0) to have infinite magnitude or to 
be indeterminant. For such u{0): 

• If guessDE (u, z) > 0, then we can request o [z^), truncating if the resulting domi- 
nant exponent is actually negative. 

• Otherwise we can request an iterative determination of the appropriate order, with 
the proviso that the order be o (z'^) if a is actually positive, or else with c not a 
branch point. 

For a product of two or more operands, let u contain a proper subset of the operands 
and V contain the complementary subset. Use recursion if u and/or v is thereby a 
product. Letting U and V denote the corresponding truncated series, we can guess 
initial dominant exponents a and /3, then alternatively increase them if necessary until 
either U or V reveals its true dominant exponent. Then we know enough to compute 
the other operand series to the appropriate order without iteration, after which we know 
enough to truncate or compute additional terms of its companion if necessary. We can 
terminate and return for both U and V^ if a guess for f3 yields U = O + o (z^~^) , a guess 
for a yields U = + o (2'^"°) , and a + (3 > k. Algorithm 1 presents details. 

The chances for needing to truncate or iterate are reduced if we choose for v a factor 
for which the guess for /3 is most likely to be accurate, such as a linear combination of 
powers of z. To aid this choice we can have guessDE (...) also return a status that is an 
element of 

{lowerBound, exact, upper Bound, uncertain] , 

with these constants being a guarantee about the guess. Example rules for computing 



and propagating such a status are: 



statusOfGDE {z, z) 
statusOfGDE (m + t>, z) 

StatusOfGDE (u^ z) 






equal, 

lowerBound, 
rstatus0fGDE(u,2) 

upperBound 

lowerBound 



iik>l, 

if StatusOfGDE (u, z) 

if StatusOfGDE (m, z) 



lowerBound, 
upperBound, 



statusOfGDE iu, z) otherwise. 



For example, if statusAndGDE {u, z) returned [lowerBound, 6] and statusAndGDE (t>) 
returned [exact, 4], then we know that series {uv, z = 0,o (z^)) is + o [z^) without com- 
puting series for u and v. 

We can return and exploit a similar status for guessinc (...). For example, if 



statusAndGDE {u, z) — > 
statusAndGlnc iu, z] — > 



[lowerBound, 6] 
[exact, 0] , 



then there is no point to iteratively increasing the requested order for u beyond 6. 

Another way to compute operand series to the necessary and sufficient order, pio- 
neered by Norman |8], is to use lazy evaluation, streams, or Lisp continuations. The idea 
is to generate at run time a network of co-routines that recursively request additional 
order or additional non-zero terms for sub-expressions on an incremental as-needed basis. 
The above guesses and iterative techniques are relevant there too, because if the request 
is for an increment to the order, it might not produce another non-zero term and if the 
request is for an additional non-zero term, iteration might be necessary to produce it. 
However, for such algorithms that don't recompute all of the terms each iteration, it is 
probably more efficient not to increase the increment each iteration. 

2.3 Optional requested number of non-zero terms 

Often rather than a requested order, users need a requested number of non-zero terms, 
regardless of what order is required to achieve that. Most often the needed number 
of terms is 1. For example, a particularly effective way to compute many limits is to 
compute the limit of the dominant term. This is particularly helpful for indeterminacies 
of the form oo — oo. As another example, if we equate a truncated series to a constant, 
then it is much more likely that we can solve this equation for z if there are only one or 
two terms in the truncated series. Thus it is important to implement a separate function 
such as 

nTerms {expression, variable = point, uumherOfNonZeroTerms) . 

Parameter uumherOfN onZeroTerms could default to 1 and/or there could be a sepa- 
rate function such as 



dominant Term {expression, variable = point) 



For a hierarchical series the user often doesn't know a priori an appropriate set of basis 
functions for the series, and the dominant basis function can be an essential singularity 
or a logarithm rather than z. For such series it is much more appropriate for users to 
request the desired number of terms rather than an order in z. In this context, it is 
appropriate to count recursively-displayed terms of a hierarchical series as if they were 
fully expanded. For example, only one such distributed term is necessary for purposes 
such as computing a limit. 

It is dangerously misleading to include a term unless all of the preceding terms are 
fully developed, which might and often does require infinite series for the coefEcients of 
some preceding terms. For example, it is inappropriate to include the z^ term of 



--+ E^h' + ^-'H-"!-') w 




if the series for the coefficient of z"^ is truncated, because (In 2;) z'^/2'^^^ dominates z^ 
for all /c > 0. This is another reason that a requested distributed term count is more 
appropriate than an order request for truncated hierarchical series. This is also a good 
reason for providing the option of not expanding coefficients as sub-series where the 
implementation can't express them as infinite series and they don't terminate at a finite 
number of terms. For example, there should be an option for even a hierarchical series 
function to return 

for expression (jl]) 

To achieve a requested number of non-zero terms, we can iteratively increase the 
requested order until we obtain at least that number of non-zero terms, then truncate 
any excess terms. The iteration could begin with a requested order somewhere in the 
interval guessDE (m, 2;) -|- [0, (n — l)guessInc(M, 2;)], where n is the requested number of 
non-zero terms. If this attempt exposes no terms, then we can increment the request by 
n ■ A where A starts at guessInc(M, z) and doubles after each failed attempt. 

However, the implementation should address the fact that an expansion might ter- 
minate as exact with fewer than the requested number of non-zero terms. The fact that 
the returned number of terms is less than requested is a subtle indication that the series 
is exact, but an explicit error-order term of the form ((. . .)°°) makes this fact more no- 
ticeable. This is additional motivation for having each intermediate series result include 
an indication of exactness, if known, as elaborated in subsection 13.21 

3 Issues about displaying truncated series results 

In contrast to most other expressions, the terms of a truncated Puiseux series expanded 
about w = Wo for finite wq are traditionally displayed in order of increasing powers of 
w — Wq, even if there are logarithms involving w in expressions multiplying some of those 
powers. If series are represented using the same data structures as general expressions 
but different ordering, then the different ordering might prevent key cancellations because 

10 



efficient bottom-up simplification typically relies on the simplified operands of every 
operator having the same canonical ordering. For example, 



~z + series (e^, z = 0, o (-2^)) — > — 2; + I 1 + 2; 



z" 



-^ -z + l + z + ^ 



2 

2 



rather than simplifying all the way to 1 + z^/2. 

The series (. . .) function is most often used alone as an input, perhaps with the result 
assigned to a variable, rather than embedded in an expression. If so and the assigned 
value is ordered normally or the system re-simplifies pasted and assigned values when 
used in subsequent expressions, then the differently-ordered series result would safely be 
re-ordered into the non-series order during simplification of that subsequent input. 

One way to overcome this difficulty entirely is to use a special data structure for 
series results, then use a special method for displaying those results. However, the next 
subsection describes how onerous it is to integrate such special data structures into a 
system in a thorough seamless way. 

3.1 The pros and cons of an explicit error order term. 

Maxima 5.18.1 displays "-I-..." at the end of a truncated series result, which means 
o (z") where n is the order argument provided by the user. Mathematica 7.0.1.0 displays 
a big-0 term, and Maple 13.0 displays a big-0 term if the result isn't exact. Even when 
result orders are always as requested, a displayed ellipsis is useful and a displayed error 
term is even more valuable. They remind users that although the result is symbolic, 
it is perhaps or definitely approximate. Moreover, it provides an opportunity for the 
implementation to make such truncated series infectious, which helps prevent users from 
misusing inappropriate mixtures of approximate results with exact results or with results 
having different orders or expansion points. For example, if a result of series (...) is 

{w - 27r)"'/' + iw- 2nf + o {{w - 27if) , 

then adding sin (w) to this result would return 

[w - l-nY^I^ + (w - 27r) ^[w- luf + o ((w - 27r)^) . 

With this infectiousness, / {yS) + o [{w — wq)™) is an elegant and convenient alternative 
to the input series (/ [w) , w= Wq, o {{w — Wq)"^)). 

Unfortunately, the effort required to do a thorough job of implementing this syntactic 
sugar is extensive. To correctly propagate the influence of an error order term, every 
command, operator and function should have a method for properly treating it. This 
obligation also applies to every new command, operator or function that is subsequently 
added to a system, including user-contributed ones that aspire to be ffist-class citizens 
seamlessly integrated into the system. 

Table |2] shows some examples that test an implementation's handling of explicit error- 
order terms. Also, test if the implementations you use can directly plot series results or 

11 



apply operators and functions such as J, ^, lini, solve (...), and series (.. .) to series 
results without the nuisance of first explicitly converting the series result to an ordinary 
expression. 

3.2 Computation, propagation and display of an order term 

With a little-o interpretation of the series (...) function order parameter and strict ad- 
herence to delivering as-requested order, it is unnecessary to represent and propagate 
error-order during the internal calculations, even if we display, o (2;^) for that requested 
order. 

Mathematica 7.0.1.0 and Maple 13.0 display an error order using O rather than o. 
However, correctly determining a correct and satisfying exponent to use in O requires 
more work than o: To return a result with a requested order 0(2;^) using a O {z'^), we 
must determine a z/ > r such that the degree of the first omitted non-zero term, if any, 
is at least z/. We can't just display O {z'^^^), because with fractional powers the exponent 
of the first omitted term can be arbitrarily close to r. 

One way to determine a correct r is to actually compute the first omitted non-zero 
term, but not display it. However, that omitted term can have a degree arbitrarily 
greater than r, costing substantial computation. Moreover, there might not be any non- 
zero terms having degree greater than r. Therefore we don't know in advance what order 
if any will just expose a next non-zero term whose coefficient we discard. Also, if we find 
such a term, it would be more informative to display O {z'^) rather than O {z^)- 

Another way to determine a correct z/ is to compute a series that is o (2;'^+^) with a 
predetermined A such as 1, then truncate to o{z'^) and display O {z'^) or {z'^) where 
u is the dominant degree of the truncated terms, or display O [z^~^^) if no terms were 
truncated. However, with fractional exponents there can be arbitrarily many non-zero 
terms having exponents in the interval [t^,t -\- A], which is costly. Moreover, users 
might judge the implementation unfavorably if the exponent in O is obviously less than 
it could be. For example with sin z the series z — z^/3 + O (z^) is disturbing compared 
to z — z^/3 + O (z^), which can be more informatively displayed as z — z^/3 + O (z^)- 

When computing series, we often know the exact order for the series of some or all 
sub-expressions. For example, if the requested order is 3 then z'^ is (-2°^), whereas 
z^ is O (z^)- If we decide to store error-order information with the series for each sub- 
expression, then it preserves information to store with the error order whether it is of 
type o, O, or G, and to propagate it according to rules such as, for a < /3; 

e(z°) + 0(2^) -^ 6(2"), 
e(z") + 0(2") -^ o(z"), 
e(z") + o(2") -^ o(z"), 

Q{z") + o{z") -^ 6(2"), 
0(2") + 0(2") -^ 0(z"). 
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4 A frugal dense representation 

A sparse series representation can more generally accommodate truncated Hahn series, 
which can also have irrational real exponents. For example, 

z^ + z-^ + z^ + ... 

This extra generality is desirable for hierarchical series. Adaptive-precision interval arith- 
metic can be used to keep the exponents properly ordered. 

However, most published algorithms for series are written in a notation that encour- 
ages a dense representation as an array or list of coefficients with implied exponents. 
Adding two series is easy for sparse representation. Otherwise, adapting most of the 
published truncated power series algorithms to a sparse representation seems likely to 
make them more complicated. Moreover, with typical applications dense representation 
is efficient for most univariate polynomials, hence also for most series. As described in [H] 
recursive dense representation is also surprisingly efficient for most sparse multivariate 
polynomials, hence also for recursive hierarchical series or multi-variate series. Therefore 
this section describes a particularly efficient dense representation and some algorithmic 
necessities for maintaining it. 

The allowance of negative exponents suggests that we should also explicitly store the 
exponent of the dominant term. 

The allowance of fractional exponents suggests that we should also store the implicit 
positive rational exponent increment between successive stored coefficients. To minimize 
the number of stored coefficients, it is most efficient to make this increment be the 
greatest common divisor of all the exponent increments between successive non-zero 
coefficients. (The gcd of two reduced fractions is the gcd of their numerators divided by 
the least common multiple of their denominators). 

The truncated series can be represented canonically as a leading exponent of 0, an 
exponent increment of 0, and consistently either no coefficients or one coefficient that is 
0. Using no coefficients is more frugal and easier to program, but one zero coefficient 
permits distinguishing a floating-point series 0.0 + 0(2;'^) from a rational-coefficient series 
+ o{z^). 

Rather than the gcd of the exponent increments, many implementation instead use 
the reciprocal of the common denominator of all the exponents. However, this can require 
arbitrarily more space and time. For example, it would store and process 21 coefficients 
rather than 3 for 1 + 2z^°/^ + 3z^°/^, and it would store and process 31 coefficients rather 
than 4 for z'^^ + 2 + 3z^^. 

Either way, for canonicality, programming safety and efficiency, it is important to 
trim leading, trailing and excessive intermediate zero coefficients from intermediate and 
final results wherever practical. However, within a function that adds two series, etc., it 
might be convenient to temporarily use series that have leading coefficients and/or an 
exponent increment that is larger than necessary. For example: 

• When two series having different exponent increments are multiplied, we can use 
copies in which extra zeros are inserted between the given coefficients of one or 
both series so that their mutual exponent increment is the gcd of the two series 
exponent increments. 
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• Let 7 be the gcd of the dominant exponents and exponent increments of two se- 
ries. If the series have different dominant exponents and/or different exponent 
increments, then before the series are added, copies of one or both series can be 
padded with extra zeros before the dominant coefficient and/or between coefficients 
so that both series start with the same imphcit exponent and have the same imphcit 
exponent increment 7. 

For both examples the resuhing series should then be adjusted if necessary so that its 
leading coefficient is non-zero and its exponent increment is as large as possible. 

When computing any one series, the sub-expressions all have the same expansion 
variable and expansion point after the transformations described in Section 1. Also, 
the desired order of the result is specified by the user and can be passed into the recursive 
calls for sub-expressions, adjusted according to Table [2l If an implementation delivers 
an as-requested o-order, then there is no need to store it in the series data structure for 
intermediate series results. Therefore, only the dominant exponent, exponent increment 
and frugalized coefficient list or array are necessary for an internal data structure during 
computation of any one series. 

For each function or operator, such as In and "+", it is helpful to have a function that, 
given an expression and a requested order for expansion at 2; = 0: 

1. guesses the dominant exponents of the operand series where needed, 

2. computes the guessed necessary and sufficient order to request for the operand 
series from Table |2] and the guessed dominant exponents, 

3. recursively computes those series to the guessed necessary and sufficient orders, 

4. truncates if the requested orders are excessive, or iteratively increases the requested 
orders if they are insufficient, 

5. invokes a companion lower-level function to compute the result series from the 
resulting operand series. (For computing a function of a given series, this companion 
function would be invoked directly.) 

If we wish to report to the user the type of the resulting order {6, o, or O) and the 
corresponding exponent, then the internal data structure must also contain fields for 
those. 

If we also wish to preserve with the final result the expansion variable and expansion 
point, then we must have an external data structure that includes those together with 
the internal data structure. 

5 Exponentials interact with logarithmic coefficients 

Definition. A function / {z) is sub-polynomial with respect to z if 

r \{ ie\^ -f- f ie\\ )^ V7 > 0, 

lim re f (re )\ = < 

r^o+ 1^ ^ ' ^ ^1 [cx) V7<0. 

Examples include 
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• an expression independent of the expansion variable z, or 

• an expression that is piecewise constant with respect to z^ or 

• an expression of the form In (c i^z) 2;"), where a G Q and c iz) is sub-polynomial, or 

• any sub-exponential function of sub-polynomial expressions. 

Definition. A sub-exponential function g i^z) is one for which g (In 2;) is sub-polynomial. 

Examples of sub-exponential functions include rational functions, fractional powers. In, 
inverse trigonometric and inverse hyperbolic functions. 

Non-constant sub-polynomial coefficients can arise for a series of an expression that 
contains logarithms, inverse trigonometric or inverse hyperbolic functions of the series 
variable. 

Most Puiseux-series algorithms require no change for coefficients that are generalized 
from constants to sub-polynomial expressions, making this powerful generalization of 
Puiseux series cost very little additional program space. For example, some algorithms 
for computing functions of constant-coefRcient Taylor series are derived via a differential 
equation. However, once the algorithms are obtained for constant coefficients, there is no 
need to incur the difficulties of including any dependent coefficients in the differentiations 
and integrations used in the derivations. 

However, the algorithm for computing the exponential of a series does require a 
change: For series having a non-negative dominant exponent, the algorithm begins by 
computing the exponential of the degree-0 term to use in the result coefficients. If the 
degree-0 term is a multiple of a logarithm of a monomial containing a power of z, then the 
exponential of the leading coefficient generates a power of z. To avoid incorrect truncation 
levels, this power of z should be combined with the (perhaps implicit) power portions of 
the data structure so that the true degree of each term is manifest in a canonical way. 

Other places where coefficients interact with powers are computing derivatives or 
integrals of a series with respect to z. For example, ^Inz — ?■ z~^, and Jlnzdz — ?■ 
{Inz — 1) z. 

6 Essential singularities 

Exponentials and sinusoids of negative powers of z are essential singularities at z = 0. For 
a series, let U_ be the sum of the terms having negative exponents and U be the sum of all 
the other terms. We could use the transformation e— """^ — )■ e—-e^, then compute the series 
for e^, then distribute the essential singularity e— over the resulting terms as factors in 
the coefficients. We could similarly use angle sum transformations for sinusoids of series 
having negative exponents. However, essential singularities dominate any power of z at 
z = 0. If we distribute the essential singularity, then subsequent series operations can 
truncate terms that dominate retained terms that don't contain the essential singularity, 
giving an incorrect result. For example, 

series fe^"' (1 + 2;^) + sin z, z = 0, 0(2;) j — > fe^"' + e''"'^;^ j + (2; + 0(2;)) 

incorrect 7^^ , , / \ 

— > e^ +z + o{z). 
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Therefore it is more appropriate to produce a recursively represented series in e— 
having coefficients that are generahzed Puiseux series in 2; - a hierarchical series. 

We could have one extra field in our data structure for a multiplicative essential singu- 
larity. However, that complicates the algorithms for very little gain, because subsequent 
operations can easily require a more general representation. For example, one field for a 
single multiplicative essential singularity can't represent 

e"" + z. 

The extra effort of implementing such a limited ability to handle essential singularities 
is better spent implementing more general hierarchical series. 

Collecting exponentials in an expression permits computation of generalized Puiseux 
series for some expressions containing essential singularities that are canceled by the 
collection. For example, 

gcot^ 2 8 

(e^/^)""" |x G M ^ e(^'°^)/^-^ei-^'/6+-^e - — + ■■■. 

6 

7 Unnecessary Restrictions 

Not all Puiseux-series implementations currently allow fractional requested order. How- 
ever, if fractional exponents are allowed in the result, then it is important to permit them 
as the requested order too. Otherwise, for example, a user will have to compute and view 
1001 terms of exp (^^1/1000^ merely to see the first two terms 1 + z^^^^*^^ . 

Not all implementations currently allow negative requested order. However, if neg- 
ative exponents are allowed in the result, then it is important to permit them as the 
requested order too. Otherwise, for example, a user will have to compute and view 1001 
terms of e^ / z^^^^ merely to see the first two terms z~^^^^ + z~^'^^ . 

These restrictions are probably caused by restricting some field in the data structure to 
a one-word signed or unsigned integer, which can also unnecessarily limit the magnitude of 
the requested order. Although most likely motivated by a desire for efficiency, the savings 
are probably a negligible percentage of the time consumed by coefficient operations and 
other tasks. 

Not all implementations currently allow non-real infinite-magnitude expansion points, 
such as for 



{w "^, w = ioo, 3) — )> w ^ + 6 {w''' 



series [w 

despite the fact that such limit points can be mapped to a real infinity by a transformation 
such as w — )■ —iz. 

Not all implementations currently allow full generality for sub-polynomial coefficients. 
For example, 

series (arcsin (In 2;) , 2; = 0, 3) — )■ arcsin (In 2;) -|- 6 (z°^) . 
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The sub-polynomial coefficient of z^ can be developed as a truncated hierarchical infinite 
series In^; + (ln2;)^/6 + ■ ■ ■ , which is preferable for some purposes such as computing a 
limit at z = 0. However, if the request is for expansion in powers of z. then arcsin (In 2;) 
has the advantage of being exact and much simpler. 

Summary 

The generalization from Taylor series to generalized Puiseux series introduces a surpris- 
ing number of difficulties that haven't been fully addressed in previous literature and 
implementations. Such issues discussed in this article are: 

• avoiding unnecessary restrictions such as prohibiting negative or fractional orders, 

• the pros and cons of displaying results with explicit infectious error terms of the 
form o (. . .), O (. . .), and/or O (. . .), 

• efficient data structures, and 

• algorithms that efficiently give users exactly the order they request. 
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Requested operand orders for a result having order o (2;'^) , with 
Table 1: f/= cz" + te'" + - ■ - + (z™) andV = az^ + hz"' + - ■ ■ + o{z'^) 

where a, (3, a, 7, m, n, /c G Q 



operation 



result dominant exponent 



request m and n 



U±V 



> min (a, /3) 



m 



n 



k 



UV 



>a + [5 



m 

n - 



■-k-l3 
k — a 



U 
V 



>a- (3 



m = k + (3 
n = k — a + 2/3 






> 7a 



m 



k + {1 — 7)0; 



e 

cosU 
coshf/ 



if a > 

-00 otherwise 



m 



k if a > 

essential singularity otherwise 



sinf/ 

tanf/ 

sinhf/ 

tanhf/ 



a if a > 
-00 otherwise 



m 



k if a > 

essential singularity otherwise 



Inf/ 



0" if cz°' = 1 
otherwise 



m 



k 



a 



arctanh U 



a if q; > 

otherwise 



m 




if cz" = 1 V C2" 
if a < 
otherwise 



arctan U 



a if q; > 

otherwise 



m 




if cz" = iy cz'^ 
if q; < 
otherwise 



—I 



arcsinh \J 



a if q; > 

otherwise 



vn 




if cz^ = i\l cz^ 
a < 
otherwise 



—I 



arcsin U 



a if q; > 
otherwise 



m 




if cz" = 1 V cz" 
a < 
otherwise 



arccos U 
arccosh U 



a/2 if C2" = 1 
otherwise 



m 




if cz" = 1 V C2" 
if a < 
otherwise 
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Algorithm 1 Compute series(-u) and series(f ) for computing series(Mv) to o (z'^) 
Input: Symbolic expressions u and v, variable z, and rational number k. 
Output: The ordered pair of truncated series [U, V], each computed to the necessary 
and sufficient order so that UV is o (-2^)- 

6u < 1; S^ < 1; // -1 nieans these increments haven't yet been computed 

a ^ guessDE (-u, z) ; f3 ^ guessDE {v, z); 
m ^r- rriQ <— k — /3; n ^ uq ^ k — a; 
loop 

U ^ series (u, z = 0, o (-2™)) ; 
if f/ ^ + o (z") , then 

a -^ dominantExponent (U) ; 
n ^ k — a; 

V •(— series (f , z = 0, o {z"-)) ; 
iiV = + (2") , then return [U, V] ; 
/3 <— dominantExponent {V) ; 
ii m = k — P, then return [U, V] ; 
ii m > k — 13, then return [truncate {U, k — P) ,V]; 
return [series [u, z = 0, o (2;™"'^)) , V] ; 
V ^ series (f , z = 0, o {z"')) ; 
iiV ^0 + (2") , then 

/3 ^ dominantExponent {V) ; 
m ^ k — (3; 

U ^ series (m, z = 0, o (-2™")) ; 
if f/ = + o (2") , then return [f/, V] ; 
q; ^ dominantExponent (U) ; 
if n = A; — a, then return [U, V] ; 
if n > A; — a, then return [U, truncate {V, k — a)] ; 
return [f/, series (f, z = 0, 0(2;"^"))] ; 
if m + n > A;, then return [U, V] ; 

1 6u < A guesslnc {u, z) = 0, 

5u -^ \ guesslnc {u,z) 5„ < 0, 

5u + (^M otherwise; 

1 5t, < A guesslnc (f , 2;) = 0, 

5v -^ \ guesslnc (t>, 2;) 5^ < 0, 

5j, + (5j, otherwise; 

m ^ min (mg + 5„, A; — n) ; 
n ^ min (no + 5y,k — m) ; 
endloop; 
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Table 2: Test examples for treating explicit error order terms 



Input equivalent to 


Increasingly informative results 


z — series {z, 2; = 0, 0(2;^)) 


oiz),Oiz^),0 


In (series (e^, 2; = 0, oiz^))) 


z+o{z^),z + 0{z^),z + Q{z^) 


1 + z — series (e^, z = 0, (z)) 


oiz),Oiz^),eiz') 


series (e^, z = 0, o{z^)) — series (e^ 2; = 0, o{z'^)) 


oiz'),Oiz^),Qiz^) 


series (e^+z^, z = 0, o{z'^)) — series (e^, z = 0, o{z'^)) 


oiz'),Oiz^),Qiz^) 


series (e^ + z"*, z = 0, o{z^)) 
series (e^, z = 0, {z'^)) 


l + o{z^), 1 + 0(^3), 1 + 0(^3) 
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